Small-scale structure of nonlinearly interacting species advected by chaotic flows 
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We study the spatial patterns formed by interacting biological populations or reacting chemicals under the 
influence of chaotic flows. Multiple species and nonlinear interactions are explicitly considered, as well as cases 
of smooth and nonsmooth forcing sources. The small-scale structure can be obtained in terms of characteristic 
Lyapunov exponents of the flow and of the chemical dynamics. Different kinds of morphological transitions are 
identified. Numerical results from a three-component plankton dynamics model support the theory, and they 
^— ■») . serve also to illustrate the influence of asymmetric couplings. 
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] The transport of reacting substances by a fluid flow is a problem appearing in a variety of disciplines, from classical 
^ ' combustion studies to chemical reactor design. Important environmental examples arise in the study of atmospheric ad- 
Tj- vection of reactive poUutants or chemicals, such as ozone, or in the dynamics of plankton populations in ocean currents. 

The inhomogeneous nature of the resulting spatial distributions was recognized some time ago, but more recently satel- 
^ Ute remote sensing and detailed numerical simulations identify filaments, irregular patches, sharp gradients, and other 

Q complex structures involving a wide range of spatial scales in the concentration patterns. We analyze here the smaU-scale 

structure of a large class of models of transported reacting substances in terms of basic concepts from dynamical systems 



theory, and apply the results to a model of plankton dynamics. 



, I. INTRODUCTION 

cn 

Tracers stirred by fluid motion are known to develop strong inhomogeneities, usually in the form of filamental features, arising 
from a kind of variance cascade from the forcing scale towards smaller scales These structures are observed at size ranges 
^— I as diverse as the ones relevant for laboratory experiments atmospheric transport [^|^, or temperature or chlorophyll 

patchiness in the ocean They provide sensitive mixing mechanisms and they are in some sense "catalists" for chemical 

i or biological activity occurring in the flows [p|jlo|]. For example, in the case of atmospheric chemistry, the presence of strong 

• ^ ' concentration gradients has been shown to have profound impact on global chemical time-scales [[Till. The same phenomenon 
^ has been observed in models of ocean plankton dynamics [[l2|]. On the other hand, chemical reactions or biological competition 

occurring in the advected substances modifies the characteristics of its spatial distribution [[l3|-p^. Thus the interplay between 

• ^ fluid flow and its chemical activity is an important issue to understand the spatial structure of concentration fields both in 

' environmental and in artificial flows. 
^ . For the case of inert transported tracers (also called the passive scalar problem) much progress has been achieved in the last 
■ ■ ■ ' years in describing, at least in statistical terms, its spatial structure [|l6[jl^]. Most of the results have been formulated in terms of 
structure functions, which describe the statistics of spatial fluctuations, or of variance power spectra. A regime of fluid motion for 
which considerable understanding has been attained is the so called Batchelor regime [ [l8| ] , or viscous-convective range, which is 
the range of scales in a turbulent flow below the Kolmogorov scale for the velocity field (the flow is thus dominated by viscosity 
and spatially smooth) but above the scale for which diffusion effects smooth out the structure of the transported substance. This 
regime is of appreciable extent for large Schmidt number (Prandtl number if the transported substance is temperature). It is by 
now clear that transport behavior in this range is equivalent to the one generically obtained under the conditions of Lagrangian 
chaos or chaotic advection [|T9[-pl]], that is, stretching and folding of fluid elements following chaotic trajectories in laminar 
flows. Simple two-dimensional flows, and even steady flows in three dimensions, lead generally to this kind of fluid particle 
motion. In addition to modeling laboratory situations, the chaotic advection paradigm has been shown to be a useful approach to 
understand geophysical transport processes at large scales [^]. In the Batchelor or chaotic advection regime, in closed flows, the 
power spectrum presents a k^^ decay at large wavenumbers 1 18 2^, and the structure functions behave logarithmically 25 1. 



In open flows the singularities are restricted to fractal sets of zero measure [ |2q , [27[ ], but they can be experimentally visualized | 
and affect the global scaling behavior. 
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In the case of reacting transported substances (still passive in the sense that they do not alter the flow), simple model reactions 
such as ^ + _B — > 2B [ p8| ] or A + B — > C were studied in both closed and open [ p^ flow situations. Possibly, one of the simplest 
chemical reaction models one can consider is the first-order reaction, or linearly decaying tracer. It consists simply in the decrease 
of the concentration of a chemical at a rate proportional to the same concentration. This dynamics describes the consumption of a 
reactant in a binary reaction, when the other reactant concentration is kept constant, or the spontaneous decay of a finite-lifetime 
substance, such as a radioactive tracer or a fluorescent dye. In the presence of a source the concentration reaches a non-zero 
statistically steady state. Sea-surface temperature relaxation towards atmospheric values [p9|], and the relaxation of plankton and 
nutrient concentrations towards mixed-layer values [ p^ have also been modeled in this way. Ref. [ pO[ ] pointed out the relevance 
of this model for vorticity dynamics in two-dimensional turbulence in the presence of drag. Recent results p^ , ^3pl] , p4[ ] have 
completed the classical ones by Corrsin [ ^5| ] so that we have now a rather complete picture of the spatial structure of finite- 
lifetime substances transported by chaotic flows. The decay-time constant makes the tracer power spectrum steeper than the 
Batchelor law obtained for the passive tracer, and in general the scaling of the structure functions is controlled by the relationship 
between the decay-time constant and the Lagrangian Lyapunov exponent of the flow. The basic physical mechanisms behind 
these results are the compression of fluid elements along contracting directions of the flow, with the consequent increase in the 
gradients, and the tendency of the decaying chemistry to relax these gradients. The competition between these processes leads 
to the appearance of morphological transitions as the value of the two time scales changes: when the chemical decay is faster 
than the compression by the flow, the concentration pattern is differentiable and characterized by a smooth structure. In the 
opposite case the concentration develops a rough (non-differentiable) structure which reflects in non trivial scaling exponents for 
the structure functions. The rough distributions have filamental aspect, since there is always an expanding direction for the flow 
along which the concentration is smoothed out. Thus these transitions have been termed smooth-filamental transitions [^. 

This simple picture should be completed with the recognition that stretching (and compression) is usually not homogeneous, 
and different parts of the flow experience different stretching histories characterized by the probability distribution of the La- 
grangian finite-time Lyapunov exponents of the flow. This has as a consequence that the set of structure functions display 
anomalous (multifractal) behavior Different points in the system may display different scaling behavior, but the mor- 

phological transitions mentioned above can be still identified as the change of behavior in a macroscopic portion of the system, 
that is, a part with fractal dimension equal to the full spatial dimension. The anomalous scaling is particularly pronounced in the 
case of open flows with a localized mixing region. 

The presence of the above intermittency corrections to simple scaling behavior, both in inert and in reacting flows, has a 
particular status: on the one hand it implies that some of the assumptions implicit in the simplest theoretical approaches are 
incorrect at a fundamental level; on the other hand, the corrections in quantities most directly accessible to experiments, such 
as low-order structure functions, could be rather small. Thus, it is usually a rather good approximation to neglect, in a first 
approach to the problem, multifractal behavior and concentrate into the bulk, dominant behavior. This is specially pertinent in 
environmental flows, where precise quantitative measurements are not always possible or reproducible. Comparison between 
structure functions of different orders has been however performed, and multifractality clearly demonstrated there [ p^ , ^7[ ] . 

In addition to its intrinsic value as a model for several relevant chemical situations, the linearly decaying model has been 
pointed out to represent a much broader class of chemical reactions [ ^ , ^8| |: essentially any chemistry (or biology) leading 
to a negative Lagrangian chemical Lyapunov exponent, as long as we consider just the small-scale structure. The Lagrangian 
chemical Lyapunov exponent is defined as the average rate of convergence or divergence between concentration values present 
in a particular fluid particle, when it is initialized with slightly different initial concentration values. It turns out that it is 
enough to replace the decay rate of the linearly decaying chemical model by the largest (less negative) of the chemical Lyapunov 
exponents to obtain the small-scale characterization of the structure of a nonlinear chemical or biological model. Although 
some qualitative numerical checks of this equivalence have been already presented [ ^9| ] both for open as for closed flows, no 
quantitative evaluation of the predictions of the theory has been presented so far 

In this Paper, we consider the scaling behavior of reacting fields advected by Lagrangian chaotic flows. We consider mul- 
ticomponent nonlinear chemical or biological models, and generalize previous results to include the possibility of non-smooth 
source terms. The concept of smooth-filamental transition is revisited in the presence of these non-smooth forcings, which lead 
to different kinds of patterns. The theory is tested with numerical simulations of a model of oceanic plankton dynamics in a 
simple two-dimensional closed flow. Part of our interest in addressing multicomponent nonlinearly interacting populations arises 
from the need to clarify differences which seem to be observed in the scaling behavior of different plankton species stirred by 
the same flow pOjp^THl^l ] . Our results indicate that such differences may arise, but, as far as only the small-scale behavior and 
chaotic (non-turbulent) advection are considered, they require rather asymmetrical couplings. In the theoretical developments 
we address the main scaling behavior and neglect any intermittency correction. This allows us to concentrate into the different 
morphological transitions, which we believe are the most relevant predictions to be compared with experiment or observations. 
The comparison with the numerical results is performed however in terms of the first-order structure function. The agreement 
between the theory and the numerics confirms that multifractal corrections to this low-order structure function are weak. 

The Paper is organized as follows: After this Introduction, we present in Sect. p| the general results for the smaU-scale 
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structure of the advected fields, obtained from consideration of the Lagrangian properties along particle trajectories. Particular 
plankton and flow models are presented in Sect. which also contains numerical results which support the theory and illustrate 
the variety of transitions arising from the interplay between reaction, advection and forcing. We conclude in Sect. IV with a 
summary and open questions. 



II. THE SMALL-SCALE STRUCTURE OF ADVECTED CHEMICALS 

A. Evolution models and their Lagrangian representation 

The spatiotemporal evolution of reacting fields is determined by advection-reaction-diffusion equations. Advection because 
they are under the influence of a flow, reaction because we consider species interacting with themselves and/or with the car- 
rying medium. Diffusion because turbulent or molecular random motion smoothes out the smallest scales. For the case of an 
incompressible velocity field v(r, <), the standard form of these equations is 

+ v(x, t) ■ Vc(x, t) = F(c, x) + DVM^, t) . (1) 

We consider our system to be confined in a d-dimensional box. The velocity field is described by the rf-dimensional vector 
field v(x, t). The particular numerical example of Sect. |l| will be a time-dependent smooth two-dimensional flow. c(x, t) = 
(ci (x, t ) , C2 (x, t ) , . . . , cjv (x, t ) ) is the set of N interacting chemical or biological fields advected by the flow, and {Fi, Fn ) = 
F = F(ci, ...,CAr,x) are the functions accounting for the interaction of the fields (e.g. chemical reactions or predator-prey 
interactions) including sources. The explicit dependence on the spatial coordinate x accounts for inhomogeneous distributions 
of sources and sinks of the substances, or for the spatial dependence of reaction coefficients, which act as an external forcing 
on the concentrations. If no such forcing is present, the final distributions will be finally homogeneized by diffusion, and spatial 
pattern will occur only during a transient regime. In the presence of forcing, some statistical steady state will be attained from 
the input of substances through the sources and sinks, and the transport of these inhomogeneities to smaller scales by advection. 
There is no difficulty in considering an explicit time-dependence in the sources: F = F(c, x, t). For simplicity in the notation 
we do not write explicitly this dependence, but it is easy to see that our results (in particular Eq. ( p7[ ) below) are not altered by 
this. 

The scaling properties of the concentration fields c would depend on the smoothness properties of both the velocity field v and 
of the interactions F. For the velocity field we consider a smooth spatial behavior, which corresponds to the situation of chaotic 
advection, or to the Batchelor regime in a high Schmidt number turbulent flow. In previous works [ ^2l^ ] we have assumed also 
a smooth dependence of F both on the concentrations and in space. This allowed us to write 

F(c + (5c, x + Sx) = F(c, x) + DF(c, x) ■ (5c + VF(c, x) • (5x + ... (2) 

The ellipsis denotes higher order terms. DF is the (square) matrix made of the derivatives of the N components of F with 
respect to the concentrations c = (ci, cat). VF is the (rectangular) matrix of derivatives of the components of F with 
respect to the spatial coordinates. 

There are situations in which the spatial dependence of F is not smooth. This will occur when source terms arise from 
physical processes that lead themselves to nontrivial scaling properties. For example, water temperature is an important variable 
that controls growth rates in plankton models. But it is itself a quantity which is transported by the flow and subjected to processes 
that provide to it a nonsmooth spatial structure. Then, source terms depending on temperature may be not well modeled by (^). 
Other examples have been already presented in the literature, as the case of plankton grazing on a multifractal nutrient field 
[p^]. In most cases the structure of all the fields will be determined by the combined effect of their mutual coupling and of the 
flow. But in cases in which couplings are unidirectional, that is, one of the components influences the others but there is no 
feedback on it, the influence of the component which does not receive feedback is better modeled as a source term forcing on 
the remaining components. An explicit example of this situation will be presented in Sect. 

For the case of nonsmooth source forcing, we will use 

F(c + (5c, X + (5x) = F(c, x) + DF(c, x) • (5c + S^F + ... (3) 

instead of (^. SxF = SxF{c, x, Sx) is the set of spatial increments SxF — {SxFi, SxF^), with SxFi = Fi{c, x + Sx) — 
Fi (c, x), which scale at small distances as 

\6xF,\ ^ \Sxf^ . (4) 



3 



f3i = 1 for smooth sources. In cases in which F is bounded as a function of space, as in the concrete model to be discussed 
below, the behavior (Q) should cross over a saturation behavior for increasing |(5x|. We call L the length scale at which this will 
happen (we assume it to be the same for all chemical sources labeled by different values of the index i). It should be of the order 
of the largest scale of inhomogeneity in the source terms S^F, usually of the order of the system size. The saturation value of 
\SxFi \ will be then of the order of L^^' . 

Diffusion effects are only important at the smallest scales and we will neglect them in the following. In this limit of zero 
diffusion D ^ the above description can be recast in Lagrangian form. To this end we introduce some notation: Xj^^ (xo) will 
denote the position at time i of a fluid particle that at some initial time to was at xo, that is, it is the trajectory solution of 

^Xl=v{Xl,t) (5) 

satisfying Xf°(xo) = xq. The set of Lagrangian concentrations inside this fluid element will be denoted by Cl^{xo). The 
relationships between the Eulerian and Lagrangian concentrations are 

C^xo)=c(x = X^xo),t) (6) 



c(x,t) = C*Jxo=X*°(x)) . (7) 

The last expression states that the relevant particle trajectories to recover Eulerian concentrations at position x are the ones that 
end at x at time t > to, that is, solutions of Eq. (|]) integrated backwards in time from final conditions at x. 

Equation (|l]) in the Lagrangian framework, that is the equation ruling the chemical dynamics inside the fluid element which 
foUows the trajectory X^^ (xq), reads 

|c*„(x„) = F (C*,,(xo),x = X*„(xo)) . (8) 

The set of equations (|]) and (^, which we call the flow and the chemical subsystem, respectively, are the basic starting point 
for our analysis. Note that the coupling between both equations only appears if there is space dependence in F, that is, if there 
are inhomogeneous sources. 

The quantity we are interested in is the difference in concentration between neighboring points: 

5c(x, t; 6x) EE c(x + (5x, t) ~ c(x, t) (9) 

and in particular in its scaling behavior at small distances (but larger than the diffusion scale, which we are neglecting), which 
defines the set of Holder exponents a^: 

|(5c^(x,t;5x)| ~ l^xl"' (10) 

We will consider here only the large-time statistically steady state obtained under forcing. Thus the time in ( p^ should be 
considered to be large enough for the initial concentrations at to to be forgotten. Only after this t ^ oo limit should the 
small- |(5x| behavior be considered. 

We allow for different scaling behavior in each concentration q, although we will see that this will not be the usual case. 
In general, there will be space- and time-dependent prefactors in ([l^), but one expects the exponent, as any other statistical 
characteristic of the concentration pattern, to be a constant in space and time. Under the hypothesis of our approach, this is 
indeed the case. It should be said, however, that because of the intermittency corrections we will neglect, ai may have values 
different from the constant calculated below in sets of points of zero measure. 

The Lagrangian quantity analogous to (^ is the difference in concentration between two fluid particles: 

<5C*„(xo;<5xo) = C*„(xo + <5xo) - Cj„(xo) . (11) 



B. Lyapunov characterization of the flow subsystem 

We introduce the difference between two trajectories that start at points separated by Sxo'. 

:^xo;<5xo)^X*^(xo + <5xo)-X*^ 



5X* (xo;<5xo) ^ X* (xo + <5xo) - X* (xq) . (12) 
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If this quantity is small, its dynamics is given by the linearization of Eq. that is 

|^X*„=J(XU-5X*„, (13) 

where J is the Jacobian matrix of the derivatives of v with respect to the spatial coordinates. The solution of ( pj| ) can be written 

as 

5X*^(xo;<5xo) -T(t2,ii)-<5X*J(xo;^xo) (14) 
in terms of the d x d fundamental matrix T, which is the matrix solution of 

d 

h 

with initial condition T(ii, ti) = I, the d x d identity matrix 



^^^T(t2,ti)-J(X*^).T(t2,ii) (15) 



Equation (13) is the well-known variational equation associated to the flow (g|). Some hypothesis about the flow should be 
made to obtain concrete information. A convenient assumption is to assume that (||) is a hyperbolic and ergodic dynamical 
system |p3|]. This means that at every point of the system one can identify contracting and expanding directions, associated with 
Lyapunov exponents Ai > A2 > ... > Ad, which give exponential growth or decay of |(5X(^| at long times. Common flows 
are never perfectly hyperbolic, and there are points at which, because of the presence of KAM tori, or because of tangencies 
between stable and unstable manifolds, the characteristic directions are not defined. For sufficiently cliaotic flows such situations 
are not frequent, and it is a good approximation to assume full hyperbolicity. This approximation allows considerable generality 
in the analysis, and will be implicit in the following. Since the flow is incompressible, the sum of the Lyapunov exponents is 
zero, and thus for chaotic flows, Ai > and Ad < 0. For two-dimensional flows, A2 = — I'^il- For t2 ^ ti, the singular 
value decomposition of the matrix T{t2,ti) will be dominated by the largest eigenvalues, which are related to the Lyapunov 
exponents, so that the action of T on a generic displacement vector 6x will be given by 

T{t2,ti) ■ 6x w 7ri(i2)e^^(*^~*^Vl(ti) • Jx . (16) 

TTi and ttJ are the right and left eigenvectors, respectively, associated to the singular value e^i^'^-ti) jjj physical terms, they are 
the unstable directions attached to the points x — Xj^^ and x = X*^^, respectively. Analogously, if t2 <^ ti, the singular value 
decomposition of T(i2, ii) will be dominated by the most negative Lyapunov exponent A^, and; 

T{t2,h) ■ « 7Td{t2)e^''^''-*'^7Tl{h) ■ Jx . (17) 

We note, and this will be relevant in our results, that at each time ti, Eqs. (|l6|) and ( p^ ) will be valid for all orientations of 
(5x except for orientations perpendicular to TT^^{ti). Along these directions, the action of T will be associated to subdominant 
Lyapunov exponents, as discussed below. 

If |i2 — ill — > 00, expressions ( |l6| ) or ( p^ ) would indicate that SX.^^ grows without limit. At some moment the linearized 
evolution (|l3]) will not longer be valid, and (5Xj^ will saturate at a value of the order of system size, or of some characteristic 
length scale of the velocity field. We call this length scale L. For simplicity we use the same symbol for it as for the length scale 
of saturation of S^F. If they are different, the discussion after Eq. ( P5| ) below implies that only the smallest of these length scales 
enters into the analysis. Saturation of will happen at a time t = t((5xo) such that 6^1^ (xq; (5xo) ~ L, or 

r(x„)«-llog^, (18) 

where A is either Ai > if we are using ( [T^ so that r > 0, or Ad < 0, if we are looking for the evolution towards the past (jl 
so that T < 0. 



It should be noted that both (16) and (jlTp give only the typical asymptotic behavior at large time differences. It needs to be 
corrected at least in two aspects, even within the hyperbolicity hypothesis implicit in our approach: On the one hand, at finite 
times the Lyapunov exponent has not reached completely its long-time value, but it has a value which depends on the initial 
condition and has a characteristic probability distribution [^^^. On the other hand, even at infinite times, there is a (fractal) 
set of spatial points for which the Lyapunov exponent may differ from the typical asymptotic value. This set has zero measure 
because the distribution of finite-time Lyapunov exponents becomes narrower and narrower in time, but its presence may affect 
some of the scaling behaviors described below. These two features are not independent. Both arise from the characteristic slow 
approach of the Lyapunov exponent towards its asymptotic long-time behavior [^, and their consequences are also the same: 
the introduction of intermittency corrections to the scaling behavior. Following our goal of concentrating just in bulk scaling and 
transition behavior, we do not consider in the following any correction to (|6b or (Ul). 
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C. Scaling behavior of the chemical subsystem 



From (|^, the concentration difference (jllj) satisfies 

A<5C^xo;<5xo) = 
F (C^xo + <Sxo),x = X*„(xo + 5xo)) - F (C*„(xo),x = X*„(xo)) = 

F(C*„(xo) + <5C^xo;<5xo),Xj„(xo)+<5X*Jxo;<5xo)) -F(C*,,(xo),Xj„(xo)) . (19) 
For all times during which SC^^ and S'X.l^ remain sufficiently small, we can use (|]) to get 
i 5C^xo;5xo) = 

DF(C*„(xo),X*^(xo)) •<5C*„(xo;<5xo)+4F(C*„(xo),X*^,(xo);<5X*,,(xo;<5xo)) . (20) 

This is a linear equation for (5C*^ , even though the complete dynamics (|]) may be nonlinear 

The general solution of this linear system may be written in terms of its fundamental matrix M(t2, ^i), which is the N x N- 
matrix solution of 

^M(t2,ti) = DF (CllXtl) ■ M(t2,ti) (21) 

with initial condition M(ii, ti) =1, the identity matrix. As before, the homogeneous part of the linearization (|o|), or (pT|), is 
the variational equation associated to the chemical subsystem (||). It defines a set of Lyapunov exponents, which we call the 
chemical Lyapunov exponents. They describe the sensitivity to concentration initial conditions under a fixed trajectory X^^ (xo). 
We will consider here the situation in which all of them are negative. Positive chemical Lyapunov exponents lead to strong 
divergences at small scales and in such situation neglecting diffusion may not be justified. A treatment related to the present one 
but for map models which may have positive Lyapunov exponents can be found in [^. If t2 ^ ti, the dominant term in the 
singular value decomposition of M is related to the largest (less negative) chemical Lyapunov exponent that we denote by Ac. 
The action of this matrix on a generic vector Sc of the tangent space of concentration increments would be: 

M{h,h) ■ Sc « /2(t2)e^^(*^-*i • Sc . (22) 

fi and /x^' are the right and left eigenvectors of M(t2,ii) associated to the eigenvalue e^"^^*^'"*!). As before, we neglect any 
fluctuation in the value of the chemical Lyapunov exponent. 

The quantity we are interested in is not really (pi]), but the Eulerian increments (^. The later can be obtained from the former: 

5c(x, t; ^x) = SCl (X*« (x) ; JX*° (x; <5x)) . (23) 

From (|l]) and (H): 



<5c(x,i;5x) = M(i,0) •(5co(X°(x);(5X?(x)) + / dsM(t, s) • ^F (C?(x), X?(x); (5X?(x; 5x)) . (24) 

Jo 

Here and in the following we use to — 0. Scq is the concentration difference at t = 0. The first term describes the evolution 
of the initial concentration difference under the autonomous part of the chemical dynamics, whereas the integral term describes 
the cumulative effect of the source forcing. Since we address the long-time statistically steady state, and given that M behaves 
as (^2|) with negative Ac, the initial condition will be forgotten and we need only to consider the integral term. In the absence 
of forcing, the first term in (24) would give rise to chemically decaying analogs of the 'strange eigenmodes' of l|l],^. At large t. 



the i-component of (Q) with ( 22|) reads 

5c,(x,t;5x)« rds(/x(i)),e^^(*-^V^(s)-'5.F(Q(x),X?(x);<5X?(x;(5x)) . (25) 
Jo 

- is the i-component of the vector i.e., the component associated to the chemical species c^. Generically, the scaling 
behavior of the increment ( p^ of the chemical species c;, for any i at small (5x, will be dominated by the most important 
component of S^F (the one with minimum scaling exponent, /3„i)- This behavior will be calculated in the following. There are 
however situations in which this would not be the relevant behavior: if the i-component of the vector fi{t) is vanishing, or if 
the most important component of S^F at small 6x has no projection on //.t, then subdominant terms and Lyapunov exponents 
different from Ac should be taken into consideration. Since the eigenvectors fj, and fj,^ are changing in time, this singular 
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situation will not occur generically unless some special form of the couplings in the model enforce this situation at all times. The 
particular model to be discussed in Sect. has this property. In this Section, we analyze just the generic behavior. We split the 
integral into two contributions, one during which the backwards trajectory (5Xf , and thus S^F, is growing, which corresponds to 
the time interval {t — |r((5x)|, t), and the rest of the time (0, t — \t{Sx.)\), during which the value of S^F is saturated. Since we 
are using backwards trajectories, the most negative Lyapunov exponent should be used in the expression ([l8|) for r. We now 



substitute (Q) in Eq. (|25[), and then insert the asymptotic behavior of (5X|(x; (5x) = T(s, t) ■ 6x for t ^ s (Eq. (|17[)). The result 
is 

Jo ■/t-|T((5x)| 

We have omitted a number of space- and time-dependent factors. They give space and time dependence to the prefactors present 
in 5ci, but they do not affect the scaling behavior. To analyze the small scale behavior, we first take t ~> oo (changing variables 
to It = i — s is useful for this) and then analyze the scaling behavior of each term for small \Sx.\ : The first integral always behaves 
as |(5x|l'^^l/l'^''l whereas the second one has also this behavior if |Ac| < \^d\Pm, and j^xj'^" otherwise. The final result for the 
Holder exponents ( [lo| ) is then 

ai = min(/3™, |Ac|/|Ad|) (27) 



We recall here that /3m — niin(/3i, /32, (3n)- For differentiate sources (f3m = 1) we recover the previous result [32,33|]. It 
is also the result for a linearly decaying chemical, if |Ac| plays the role of the decay rate. An explicit time-dependence in the 
source terms would affect the value of Ac, but expression (^tJ) remains valid. 

As stated before, the scaUng exponent ( ^ ) will be the one obtained for generic orientation of the displacement Sx.. But at 
each point x there will be orientations for which the subdominant Lyapunov exponent A^-i should be used instead of A^: the 
directions orthogonal to the local most contracting direction tt^^. Again, in directions perpendicular to both 7r|j and tt^ j^, Xd-2 
should be used, and so on. It is enough to replace jAdl by — A in (p7[), where A is the relevant Lyapunov exponent, to get the 
scaling behavior along these directions. Along the expanding directions, i.e. the directions for which the pertinent Lyapunov 
exponent is positive, we have = /3m, which for smooth sources mean smooth scaling behavior. We note however that source 
terms having particular anisotropic properties require special consideration. This is the case of the model in Sect. |l|. 

In the same way, at each point there are particular directions 6c in concentration space such that a subdominant chemical 
Lyapunov exponent should be used instead of Ac. But these dkections will generically not be aligned with the concentration 
coordinate axes, and thus would not be relevant to the scaling of the real chemical species Sci, unless some particular form of the 
coupling is present. For the generic case, ( p7| ) states that the small-scale Holder exponent of all the interacting chemical species 
are the same. 

The minimum condition in the equations for the Holder exponents give rise to interesting morphological transitions as the 
model parameters are varied, as one or the other of the expressions under the minimum function become the relevant one. 
Examples of the transitions are given in the next Section. 



III. NUMERICAL RESULTS FOR A PLANKTON MODEL 



In the numerical investigations below we will consider a simple model of plankton dynamics stirred by a two-dimensional 
time-dependent flow. 

This plankton model is a typical predator-prey system where three trophic levels are considered: the nutrients, 
parametrized by the carrying capacity C of the water parcel (defined as the maximum phytoplankton content it can support 
in the absence of grazing), the phytoplankton P and the zooplankton Z biomass concentrations. The dynamics of these species 
is given by 

_ = Fc = a(Co(x)-C) (28) 
'^^ =Fp=p(l-^) -PZ (29) 



dt V c*, 

^ = Fz=PZ-bZ^. (30) 
at 

The Lagrangian 'chemical' subsystem (||) is simply obtained by considering that all the reactions ( psfpol ) occur in a particular 

4o 



fluid element, so that the spatial coordinate x in (|28|) is taken to be the position of this fluid particle in the ocean flow. 
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All terms have been adimensionalized to keep a minimal number of parameters. The first equation (|8]) states that the carrying 
capacity adapts (at a rate a) to the local value of some source of nutrients Cq . This will be the only explicitly inhomogeneous term 
in the model, and describes a spatially dependent nutrient input arising from some oceanic topography-determined upwelling 
distribution or latitude dependent illumination, for example. The first terms in Eq. (|9|) describe phytoplankton logistic growth, 
whereas the last one models predation (grazing) by zooplankton. This effect gives also rise to the first term in (|30|). The term 
containing b, the zooplankton mortality, describes zooplankton death due to higher trophic levels. For the parameter values we 
are using, in the absence of stirring by the flow, the system evolves to a stable equilibrium state which is non-uniform in space 
because of the inhomogeneous source Co(x). We use for the source the expression Co{x, y) — 1.7 + 0.52 sin {2tix) cos (27rj/). 
This is a smooth function and then /3c = 1- 

This particular plankton dynamics has not been chosen because of some particular biological relevance, but because it allows 
us to illustrate a rather rich variety of transitions in the context of the theory of Sect. |^ In addition, it contains some of the 
nongeneric features that are difficult to discuss in general, but that can be easily addressed in each particular case by simple 
extensions of the theory of Sect. ||, and that may be relevant for the understanding of real flows. We note that the dynamics 
of the carrying capacity Eq. ( ^8| ) influences other species, but no feedback from P or Z affects the dynamics of C. This 
leads to a linearized dynamics described by a matrix DF which has a box structure for all times, and thus, there is always 
a Lyapunov exponent, of value exactly —a, associated to a contracting direction along the direction of the carrying capacity 
fi = fj,^ = (1,0,0), which is decoupled from the ones associated to the P-Z remaining subsystem. This was one of the 



nongeneric situations discussed in Section QC. We see that it can arise rather easily in the presence of sufficiently asymmetric 
couplings between the variables. 

One way to analyze our model is to consider first the dynamics of the carrying capacity. It is independent of the other variables 
and its small-scale structure, characterized by the Holder exponent ac, will be given by the interplay between the flow Lyapunov 
exponent, the decay rate a, and the source smoothness degree Pc ~ 1> with the result (p7[): 

The influence of C into the remaining P-Z subsystem appears only through the denominator in (|9|). We can consider the 
P-Z subsystem as a 2 x 2 chemical dynamics forced by a source term C(x, t). The associated components of S^F behave as 
\SxFp\ « \SC\ w |(5x|°^ and jf^ajF^j = 0. We thus have that the smoothness exponent of the forcing into this subsystem is 
/3^^ — ac- From (P7[), and since the variables P and Z are coupled in a generic way, the Holder exponents describing the 
small-scale structure of P and Z are equal and given by 

apz = min ( ac, tt^] ■ (32) 



|A.| 

Ac is now the largest Lyapunov exponent associated to the forced P-Z subsystem. As A^, it can only be estimated numerically. 

An alternative way to analyze our model is to recognize that the theoretical arguments of Sect. || imply that ( p7| ) is valid if the 
minimum condition for every i is taken over the set of Lyapunov exponents which affect to this particular species Ci . Considering 
thus the three-component system globally forced by the smooth source Co, the same results are obtained. 

We define our system to be the unit square with periodic boundary conditions, we use the following model flow: 

Vxix,y,t) = — i^-t mod Tj cos(27r?/) 



2U ( T 

Vy{x, y, t) = — ( ^ mod ^ ~ ) cos(27ra;) . (33) 



Q{x) is the Heavyside step function. Note that the largest lengths of both the source and the flow are given by the system size, 
L = 1. The flow is time-periodic with period T. In our simulations we use U — 1, which produces a single connected chaotic 
region in the advection dynamics. Full hyperbolicity is not garanteed, but the flow is chaotic enough to apply the theory of Sect. 
H It is easy to show that the Lyapunov exponent is inversely proportional to T. Since we are in two dimensions, Ai — | A2 1. The 
numerically determined value is Ai — IA2I « 2.35/T. We also fix 6 = 0.05, and vary the values of T and a. 

A quantity that is usually considered in scaling studies of the properties of fluid patterns is the n-order structure function. It is 
defined as 

Sl{Sx) = {\6c,{^,t;S^)n , (34) 

were the average is over space. In our numerical simulations, we perform the average in ( ^4) ) over the points in 50 rectilinear 
segments across the system, or 'transects'. 
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In general, the behavior of S'* for small Sx would be of the form 

SliSx) - \6xf^ , (35) 

If all the points of the system have the same value of the Holder exponent for the chemical species i, then it is clear that ~ qai. 
In general, intermittency will introduce corrections so that will be a nonlinear function of q. For the structure functions of 
lower order these corrections are usually small [|3]], so that for example for the first order structure function we expect 

Sl{5x) ^ \6x\''^ (36) 

to be a good approximation, where the are given by ( pH ) and (^2|). 

In order to compare theory and numerics, one needs to calculate the chemical Lyapunov exponent associated to the subsystem 
P-Z. This is done by integrating, using the same Lagrangian trajectory Xq, two copies of (p8|)-(p0|) with slightly different initial 
values of P and Z, but the same initial C. This leads to a time dependent difference between the two copies, SP and SZ, 
from which we monitor {5P^ + SZ^y^^. This is fitted to an exponential that grows or decays with a characteristic time which 
identifies Ac- In general, Ac may have an arbitrary dependence on the model parameters. It can even change sign by changing 
the characteristics of the flow an thus of the trajectory (5Xg [ p3p7| ]. But in our numerical study, for the particular flow and 
chemical subsystems, fixed values of U and b, and for all the considered range of T and a to be discussed below, the value of 
Ac was contained in the narrow interval (—0.0386, —0.0380). Thus expressions ( plj ) and ( ^2| ) become rather simple functions 
of a and T, since Ac is nearly constant, and | A2 1 is inversely proportional to T. 

For T large enough (slow flow), IA2I is small and we have that ac = apz = 1: the three concentration fields are smooth. By 
decreasing T, that is, by increasing the chaoticity of the flow, Eqs. (^Tj) and ( p^ predict two different sequences of transitions 
depending on the relationship between a and Ac, which we now describe. 

First, if a remains larger than | Ac |, the first transition encountered by decreasing T is the situation in which P and Z become 
filamental, whereas C remains still smooth. We show in Figure |l| instantaneous configurations of the three concentration fields, 
for a = 0.2, T — 20, that is, after crossing that transition [^. The predicted smooth character of C and filamental of P and 
Z is clearly observed. As expected, the features become aligned with the stable and unstable manifolds of the flow. These 
manifolds change periodically in time, following the periodicity of the flow, but the scaling properties of them do not change. 
There are differences in the details of the structure of P and Z, but Fig. ||, in which we show the first-order structure functions, 
shows that the small scale behavior is similar. The theoretical predictions for the measured value of | Ac | ~ 0.0386, ac = 1 and 
apz ~ 0.33, are also shown, in very good agreement with the numerics (least-squares fitting gives ac ~ 1.0 and apz ~ 0.35), 
despite the approximations involved. 

By further reducing T, finally C will become also filamental. Confirmation of this prediction is given in Figs. ^ and_0 for 
a — 0.2 and T = 10. The measured value of Ac is 0.0380. The configurations of C and Z are displayed in Fig. |P(the 
configuration of P is very similar to Z). The corresponding structure functions are in Fig. |^, which shows that phytoplankton 
and zooplankton have the same behavior at small scales, different from the one of the carrying capacity. The predicted exponent 
for both phyto and zooplankton is 0.16, very close to the fitting to the two numerical curves apz ~ 0.14. The prediction 
ac ~ 0.85 is more different from the observed ac ~ 0.75, but it is still close. 

Direct application of the arguments of Sect. || would indicate that the P and Z patterns are not smooth along the direction of 
the filaments, since they would inherit there the Holder exponent of the source C. This conclusion is incorrect, since the source 
is also filamental, and thus strongly anisotropic. The direction of the filaments in P and Z is the same as for C. The forcing C 
is smooth along that direction, and thus this is the behavior inherited by P and Z along the filaments. 

The second sequence of transitions one can find by decreasing T occurs when a < \Xc\- then the three fields, smooth at 
large T, should become filamental at the same value of T, and remain always with the same small-scale exponent, given by 
ac = apz = a/ 1 Ad I ■ Figure |]confirms that the transition to filamental behavior has occurred for T = 20, a = 0.025 < |Ac| ~ 
0.0385. The three structure functions, displayed in Fig. ^, have the same fitting slope « 0.20, to be compared with the prediction 
ac = apz ~ 0.21. 

IV. CONCLUSIONS AND OPEN QUESTIONS 

We have discussed the properties of the small-scale structure of forced advected chemically or biologically active substances, 
and show that they can be understood from the properties of models of linearly decaying substances, as far as the decay rate is 
replaced by the largest Lyapunov exponent of the chemical subsystem. We have extended previous results to the case of non- 
smooth forcing. The results have been applied to a three-component model of plankton dynamics, which presents a particular 
asymmetric coupling which requires special consideration. The morphological transitions predicted by the theory are observed 
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in the numerical simulations. In addition, the numerical values of the predicted scaling exponents for the first-order structure 
functions are very close to the observed ones, despite of the approximations made. This agreement would certainly deteriorate 
when considering higher-order structure functions. 

An important point to remark is that our theory only applies to the small scales of the advected patterns. Figures |[ Q and 
H do not show too strong variations in slope for different scales, but this is not guaranteed for any model. This fact should be 
taken into account when applying the results of this Paper to experimental situations and, specially, to environmental flows. 
Constructing a theory valid beyond the smallest scales remains an open challenge. 

One of the consequences of our results is that, as long as one restricts to the consideration of chaotic advection, and of the 
small-scale structure, the only mechanism leading to differences in the scaling properties of different interacting fields is the 
presence of asymmetric couplings. This could be relevant to the understanding of the differences that seem to appear in the 
scaling behavior of different plankton species in the same flow [ p4| , p^|p] ]. 

We have neglected diffusion in all our considerations. The argument was that the effect of diffusion would be a smoothing of 
the singularities below a diffusion scale, remaining unaltered the behavior above that scale. Within the same framework used 
here, this belief can be justified at least in the case of monocomponent linear models [|50[|. 

We have considered just the situation in which the largest chemical Lyapunov exponent Ac remains negative even when forced 
by stirring from the chaotic flow. The consideration of positive Lyapunov exponents remains open, specially since we expect 
diffusion to play a stronger role in the very singular configurations that can be generated. 

Finally, the consideration of multifractal behavior would be a clear improvement of the theory. It seems straightforward to 
include fluctuations in the Lyapunov exponents for the flow subsystem in the same way as in [^]. The result would be that 
the structure function exponents get a corrections which depend on the probability distribution of the finite-time Lyapunov 
exponents. It seems more difficult to include the fluctuations of the Lyapunov exponents of the chemical subsystem, since they 
are not independent variables but depend on the statistics of the flow exponents. 
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FIG. 1. Carrying capacity (a), phytoplankton (b) and zooplankton (c) configurations for a = 0.2 and T — 0. The upper panels show 
snapshots of the two-dimensional system, whereas the lower ones display the corresponding concentrations along a horizontal transect taken 
along y = 0.6. The smooth character of C, and the singular one of P and Z, is clearly seen. 
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FIG. 3. Carrying capacity (left) and zooplankton pattern (right) for a = 0.2 and T = 10. Both are now rough, but with rather different 
smoothness properties. 
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FIG. 5. Carrying capacity (left) and zooplankton pattern (right) for a = 0.02 and T = 20. Both display rough behavior of the same 
characteristics. 
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